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Abstract 

Biharmonic  equations  have  many  applications,  especially  in  fluid  and  solid  mechanics, 
but  difficult  to  solve  due  to  the  fourth  order  derivatives  in  the  differential  equation.  In  this 
paper  a  fast  second  order  accurate  algorithm  based  on  a  finite  difference  discretization 
and  a  Cartesian  grid  is  developed  for  two  dimensional  biharmonic  equations  on  irregular 
domains  with  essential  boundary  conditions.  The  irregular  domain  is  embedded  into  a 
rectangular  region  and  the  biharmonic  equation  is  decoupled  to  two  Poisson  equations. 
An  auxiliary  unknown  quantity  A u  along  the  boundary  is  introduced  so  that  fast  Poisson 
solvers  on  irregular  domains  can  be  used.  Non-trivial  numerical  example  show  the  effi¬ 
ciency  of  the  proposed  method.  The  number  of  iterations  of  the  method  is  independent  of 
the  mesh  size.  Another  key  to  the  method  is  a  new  interpolation  scheme  to  evaluate  the 
residual  of  the  Schur  complement  system.  The  new  biharmonic  solver  has  been  applied 
to  solve  the  incompressible  Stokes  flow  on  an  irregular  domain. 


1  Introduction 


In  this  paper,  we  consider  a  biharmonic  equation  defined  on  an  irregular  domain  O 
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it  is  a  bounded  open  set  in  R 2  with  a  smooth  boundary  30,  un  =  ^  is  the  normal  derivative 
of  u  on  30,  and  n  is  the  unit  normal  derivative  pointing  outward,  see  Fig.  1  for  an  illustration. 
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U  =  g2(x,y) 


Figure  1:  A  diagram  of  the  problem:  a  biharmonic  equation  on  an  irregular  domain. 


Biharmonic  equations  arise  in  many  applications.  Classical  examples  can  be  found  in 
elasticity,  fluid  mechanics,  and  many  other  areas.  In  fluid  mechanics,  the  solution  u(x ,  y ) 
of  equation  (1.1)  can  be  used  to  describe  the  stream-function  of  an  incompressible  two- 
dimensional  creeping  flow  (zero  Reynolds  number),  see  Section  5  and  [5,  15]  for  example.  In 
linear  elasticity,  u(x,y )  can  be  used  to  represent  the  airy  stress  function,  see  [32].  In  the 
theory  of  thin  plates,  equation  (1.1)  can  be  used  to  represent  a  “clamped  plate”  where  /  is 
the  external  load,  and  the  solution  u(x,y)  is  the  vertical  displacement. 

If  the  second  boundary  condition  in  (1.1)  is  replaced  by  A u\qh  =  g2{x,y),  then  the 
biharmonic  equation  can  be  decoupled  as  two  Poisson  equations  with  a  Dirichlet  boundary 
condition  on  the  same  irregular  domain 

|  A v(x,y)  =  f{x,y),  (x,y)  e  fi, 

1  v(x,y)\m  =  g2{x,y), 


A u(x,  y)  =  v(x,  y),  (x,  y)  6  fi, 

u{x,y)  1^  =  gi(x,y). 


(1.4) 


This  observation  is  one  of  basis  of  our  numerical  method  in  which  we  set  Ait]^  as  an 
intermediate  variable.  It  is  much  more  difficult  to  solve  the  problem  when  un  is  prescribed 
along  dSl. 

Various  numerical  methods  have  been  developed  for  biharmonic  equations  in  the  literature 
when  u  and  un  are  prescribed  along  dQ.  One  approach  is  to  use  a  body  fitted  mesh  with 
a  finite  element  discretization.  The  discrete  system  then  can  be  solved  using  a  multigrid 
method,  see  [1,  4,  11,  31]  for  example.  There  is  not  much  difference  between  regular  or 
irregular  domains  in  finite  element  methods  except  for  the  cost  in  the  mesh  generation  and 
the  condition  number  of  the  discrete  system  of  equations. 

There  are  limited  publications  on  finite  difference  methods  for  biharmonic  equations  on 
irregular  domains,  even  fewer  with  convincing  numerical  examples.  Most  of  the  finite  differ¬ 
ence  methods  appeared  in  the  literature  are  for  biharmonic  equations  on  regular  (rectangular 
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and  circular)  domains.  For  certain  domains,  a  conforming  mapping  can  be  used  to  solve  the 
biharmonic  equations  defined  on  the  domains  [2],  Among  a  few  finite  difference  methods  for 
biharmonic  equations  on  irregular  domains,  the  remarkable  ones  are  the  fast  algorithms  based 
on  integral  equations  and/or  the  fast  multipole  method  [8,  20,  21].  These  methods  are  most 
effective  for  homogeneous  source  term  ( f(x,y )  =  0)  and  the  boundary  conditions  (  =  9i 

and  un\ga  =  g2  in  (1.1))  are  replaced  by  the  values  of  ux  and  uy  along  the  boundary.  These 
methods  probably  still  can  be  applied  with  some  extra  efforts  for  non-homogeneous  source 
terms  and  the  essential  boundary  condition  in  (1.1).  The  implementations  of  these  methods, 
especially  when  they  are  coupled  with  the  fast  multipole  method,  however,  are  not  trivial. 

In  this  paper,  we  propose  a  finite  difference  method  for  biharmonic  equations  based  on  the 
fast  Poisson  solver  on  irregular  domains  [18],  the  embedding  technique,  and  an  augmented 
approach  for  the  decoupled  two  Poisson  equations.  The  irregular  domain  is  embedded  in 
a  rectangle  and  the  biharmonic  equation  is  augmented  with  an  intermediate  unknown  A u 
along  the  boundary  <9ft  which  is  a  one-dinrensional  quantity.  We  use  the  GMRES  iterative 
method  to  solve  the  discrete  unknown  Au|an-  Each  iteration  involves  solving  two  Poisson 
equations  on  the  same  irregular  domain,  which  we  use  the  available  package  of  the  immersed 
interface  method,  and  a  specially  designed  interpolation  scheme  to  evaluate  the  residual.  The 
irregular  boundary  is  expressed  in  terms  of  a  level  set  function.  Non-trivial  examples  show 
that  the  method  has  second  order  accuracy  in  the  infinity  norm.  The  number  of  iterations  is 
a  small  constant  independent  of  the  mesh  size.  The  proposed  method  works  for  other  type 
boundary  conditions  and  three  dimensional  biharmonic  equations  on  irregular  domains. 

1.1  A  brief  review  of  related  finite  difference  methods 

Under  the  assumption  that  u(x,y )  is  a  classical  solution  of  the  biharmonic  problems  (i.e., 
u  €  (74(ft)  P|  Cll(U)  and  u  has  piecewise  continuous  second  order  derivatives  on  <9ft).  A 
popular  technique  is  the  so  called  coupled  equation  approach, 


A  u(x,y)  =  v(x,y), 
u{x,y)  |an  =  gi(x,y), 

(x,y)  G  ft, 

(1.5) 

A  v(x,y)  =  f(x,y ), 

{x,  y)  G  ft, 

(1.6) 

v(x,y)\an  =  Au(x,y)\ 

an  -  c  (union  -g2{x,y)) , 

where  c  is  a  constant,  see  [6,  27,  28,  22]  and  others.  If  we  give  an  initial  guess  v(x,y),  then 
an  iteration  can  be  generated  until  u(x,y )  converges.  In  [22],  the  coupling  constant  c  is 
taken  such  that  0  <  c  <  2zq ,  where  zq  is  the  smallest  eigenvalue  of  the  Dirichlet  eigenvalue 
problem.  When  ft  is  a  rectangular  domain,  this  formulation  can  lead  to  an  iterative  scheme 
which  converges  for  all  sufficiently  small  values  of  c.  In  this  case,  the  two  Poisson  equations 
can  be  solved  by  a  fast  Poisson  solver,  for  example,  the  Fishpack  on  rectangular  domains  [30]. 
However,  for  an  irregular  domain,  it  is  also  challenging  to  solve  Poisson  equations  efficiently. 

One  can  also  try  to  discretize  (1.1)  on  a  uniform  grid  directly.  The  classical  13-point 
stencil  for  the  biharmonic  operator  is  most  easily  derived  by  applying  the  standard  5-point 
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Laplace’s  operator  twice 


9  ^ 

— L^L^Uij)  —  i  y  4-  Ui—ij  4-  Wjj+i 

4~  v-ij—i)  4"  2(uj_|_iJ+i  4~  rtj— 14+1  “I-  lj— 1  8“  ^j+ij— 1) 
+  (ui+2,j  4 -  Ui-2,j  4-  Uij+2  4-  — 2)) - 


(1.7) 


The  local  truncation  error  is  order  of  h2.  The  finite  difference  approximation  above  needs 
to  be  modified  at  grid  points  near  the  boundary.  One  popular  choices  is  the  quadratic 
extrapolation  where  the  normal  derivative  boundary  condition  at  the  grid  points  near  the 
boundary  are  used  to  extrapolated  the  ‘missing’  (exterior)  point  in  the  13  point  stencil.  This 
results  in  a  stencil  of  the  form: 


A  Uij  — (2\ui  j  8(iq_|_i  j  4-  Ui—ij  4*  1 

4-  Uij— 1)  4-  2('Uj_|_ij+i  +  Ui—  lj+i  4-  Ui-ij-i  4-  Wj+ij— 1) 
4-  («j+2j  4-  Uij+ 2  +  Uij- 2))  +  2hun(xi-i,yj), 


(1.8) 


when  the  irregular  point  Ujj  is  adjacent  to  the  left  boundary. 

Glowinski  and  Pironneau  [7]  made  the  observation  that  the  13-point  finite  difference 
scheme  combined  with  a  quadratic  extrapolation  formula  near  the  boundary  is  equivalent  to 
solving  the  biharmonic  equation  using  a  mixed  finite  element  method  with  piecewise  linear 
elements. 

Many  similar  modifications  are  discussed  in  [10].  Proper  treatment  of  the  points  near 
the  boundary  remains  a  challenging  problem  with  these  schemes  since  inaccurate  boundary 
approximations  may  affect  the  accuracy,  but  too  complicated  boundary  approximations  may 
destroy  the  matrix  structure. 

There  are  other  alternative  finite  difference  approximations  for  the  biharmonic  operator. 
Certain  second  and  fourth-order  finite  difference  approximations  for  the  biharmonic  equation 
(1.1)  on  a  9-point  compact  stencil  are  given  in  [29].  The  approach  there  involves  discretizing 
the  biharmonic  equation  (1.1)  using  not  only  the  grid  values  of  the  unknown  solution  u(x,y ) 
but  also  the  values  of  the  gradients  ux(x,y)  and  uy(x,y)  at  selected  grid  points. 

The  standard  iterative  methods  suffer  from  slow  convergence  when  used  to  solve  the 
system  of  finite  difference  equations  for  biharmonic  equations,  see  for  example,  [9].  Direct 
solvers  can  only  be  used  for  relatively  coarse  grids.  Bjprstad  [24]  had  introduced  a  new 
iterative  method  that  requires  only  0(N2)  arithmetic  operations  on  a  rectangular  domain, 
where  N  is  the  number  of  grid  lines  in  one  coordinate  direction. 

Due  to  the  difficulties  in  handling  the  finite  difference  approximation  close  to  a  curved 
boundary,  all  the  finite  difference  methods  discussed  above  are  restricted  to  rectangular  do¬ 
mains.  It  is  the  purpose  of  this  paper  to  provide  an  efficient  finite  difference  method  for 
biharmonic  equations  on  irregular  domains.  An  advantage  using  a  Cartesian  grid  instead  of 
a  body  fitted  grid  is  that  there  is  almost  no  cost  in  the  grid  generation  which  is  significant  for 
free  boundary  and  moving  interface  problems  that  involving  solving  a  biharmonic  equation. 

The  paper  is  organized  as  follows.  In  Section  2,  we  introduce  the  main  algorithm  and 
the  fast  Poisson  solver  for  irregular  domains.  The  interpolation  scheme  to  evaluate  the 
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residual  of  the  Schur  complement  system  is  explained  in  Section  3.  Numerical  examples  with 
grid  refinement  and  efficiency  analysis  are  presented  in  Section  4.  The  application  to  the 
incompressible  Stokes  flow  is  presented  in  Section  5. 


2  The  numerical  method  based  on  an  augmented  approach 


Consider  the  solution  of  the  following  problem 

f  A v{x,y)  =  f(x,y),  (x,y)  £  fl 

\  v(x,y)\dn  =  g(x,y), 

.  (2-9) 

J  A u(x,y)  =v(x,y),  (x,y)€Sl 

\  u(x,y)\m  =  gi(x,y). 

The  solution  apparently  is  a  functional  of  g(x,  y)  which  is  defined  only  along  the  boundary 
dQ.  We  denote  the  solution  as  ug(x,y)  to  express  the  dependency  of  the  solution  on  g(x,y). 
Let  the  solution  of  the  original  problem  (1.1)  be  u*(x,y),  and  define 

g*{x,y)  =  Au*(x,y),  (x,y)  £  d£l.  (2.10) 


Then  u*(x,y )  satisfies  the  second  Poisson  equation  in  (1.1)  with  g(x,y)  =  g*(x,y).  In  other 
words,  ug*(x,y)  =  u*(x,y)  and 


du*(x,  y) 


dn 


an 


92(x,y) 


(2.11) 


is  satisfied.  Therefore,  solving  the  original  problem  (1.1)  is  equivalent  to  finding  the  corre¬ 
sponding  g*(x,y)  and  then  ug*(x,y )  in  (2.9).  Notice  that  g*(x,y)  is  only  defined  along  the 
boundary,  so  it  is  one  dimensional  lower  than  the  solution  u(x,y). 

We  call  g(x,  y)  an  augmented  (intermediate)  variable  along  the  boundary.  The  auxiliary 
equation  then  is 


dug(x,y) 

dn 


an 


92(x,y). 


(2.12) 


Therefore  the  system  is  still  closed  since  we  have  one  more  variable  and  one  more  equation. 

The  idea  is  to  start  with  a  guess  g(k\x,  y)  as  an  approximation  to  g*(x,  y).  Once  we  know 
g(k\x,  y),  we  can  solve  the  two  Poisson  equations  in  (2.9)  to  get  the  solution  ugk\x,  y).  From 
the  residual  equation  (2.12),  we  hope  to  get  a  better  approximation  g(k+1\x,  y).  The  process 
will  become  clearer  after  we  discretize  the  system  (2.9)  and  (2.12)  in  this  section. 


2.1  The  computational  frame 

We  embed  the  irregular  domain  into  a  rectangular  domain  R  :  [a,  6]  x  [c,  d]  D  fb  The  original 
boundary  d£l  then  becomes  an  interface  within  R.  We  solve  the  two  Poisson  equations  in 
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(2.9)  on  a  uniform  Cartesian  grid 


Xj.  =  a  +  ih,  0  <  i  <  m, 

Vj  =  c  +  jh,  0  <j<n, 


(2.13) 


where,  for  simplicity,  we  assume  that  h  =  (b  —  a) /m  =  (d  —  c)/n.  The  boundary  rXl  is 
expressed  as  the  zero  level  set  of  a  two  dimensional  function  <p(x,  y )  defined  on  the  entire 
rectangular  region 


<p(x,y) 


'  >0, 

if  (x,  y)  €  R  —  ft 

=  0, 

if  (x,  y)  G  dfl, 

.  <0, 

if  (. x ,  y)  G  U. 

(2.14) 


There  are  a  number  of  advantages  using  the  level  set  method  especially  for  complicated 
geometries  and  high  dimensions.  We  refer  the  readers  to  [23,  26]  for  more  information  on  the 
level  set  method. 

The  level  set  function  is  defined  at  all  grid  points  by  ipij  =  (p(xi,yj).  It  is  important  that 
<p(x,  y )  be  a  good  approximation  to  the  signed  distance  function  at  least  in  the  neighborhood 
of  the  boundary  dCl.  If  the  level  set  function  is  not  a  good  approximation  to  the  signed 
distance  function,  then  a  re-initialization  is  needed  to  make  it  a  good  approximation  to  the 
signed  distance  function,  see  [12,  13,  14]  for  the  re-initialization  process. 

Using  the  grid  function  <pij,  all  the  grid  points  can  be  classified  as  regular  (away  from  the 
boundary  dU)  or  irregular  (close  to  or  on  the  boundary  dfl)  grid  points.  Given  a  grid  point 
( Xi,yj ),  define 


=  <Pi,j+ 1,  ifiij}, 

min  •  r  4 

Vij  —  mm{(^i— ij,  ipi+ij,  ipij-i,  ipitj+ 1, 


(2.15) 


A  grid  point  (xi,  yj)  is  irregular  for  our  problem  if 


.max,  min  /  n 

<Pij  ‘Pij  <  0 


(2.16) 


is  true.  Otherwise  the  grid  point  is  regular. 


2.2  The  orthogonal  projections  of  irregular  grid  points  on  the  boundary 

The  augmented  variable  g(x,y)  and  the  augmented  equation  (2.12)  are  only  defined  along 
the  boundary  We  need  to  discretize  the  augmented  variable  and  the  equation  at  certain 
points  along  the  boundary.  Those  points  are  chosen  as  the  orthogonal  projection  of  the 
interior  irregular  grid  points  on  the  boundary,  see  Fig.  2  for  an  illustration. 

Let  Xjj  =  ( Xi,yj )  be  an  interior  irregular  grid  point  which  means  (ptj  <  0,  and  < 

0.  The  orthogonal  projection  of  x?j  is  approximated  by  the  solution  of  the  following  quadratic 
equation: 


X*  =  x  +  a  p, 


V<£> 

iv^p 


(2.17) 
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where  a  is  determined  from  the  quadratic  equation  below: 

¥>(x)  +  |V</?|  a  +  ^  (pT  He(<p)  p)  a2  =  0,  (2.18) 

where  He(<p)  is  the  Hessian  matrix  of  the  level  set  function  <p(x,y).  All  the  quantities  are 
defined  at  the  interior  grid  point  xfj  =  ( x^yj )  (<*Pij  <  0)  and  are  evaluated  using  the  sec¬ 
ond  order  central  finite  difference  scheme.  The  orthogonal  projection  computed  using  this 
procedure  has  third  order  accuracy. 

We  will  denote  these  orthogonal  projections  of  the  interior  irregular  grid  points  by  X&  = 
(x*,y*),  k  =  1,  2,  •  •  •  ,  Nb,  and  will  omit  the  dependency  on  the  grid  points  for  simplicity  of 
the  notation.  These  orthogonal  projections  are  not  ordered  and  there  is  no  need  to  do  so. 
The  auxiliary  variable  g(x,y)  is  discretely  defined  at  X&.  The  augmented  equation  (2.12)  is 
going  to  be  discretized  at  X/,.  as  well.  We  use  the  upper  case  letters  such  as  Uij,  Vij ,  G 
for  the  discrete  approximations  at  grid  points  and  at  those  orthogonal  projections  on  the 
boundary  respectively. 

2.3  The  fast  Poisson  solver  on  irregular  domains 

The  fast  Poisson  solver  for  irregular  domains  used  to  solve  the  two  Poisson  equations  in 
(2.9)  is  based  on  the  the  fast  immersed  interface  method  (IIM)  developed  in  [17]  and  a 
modified  version  developed  in  [19,  12,  13].  The  modification  is  needed  because  the  original 
IIM  in  [16,  17]  is  designed  for  interface  problems  that  are  defined  in  the  entire  domain  with 
discontinuities  occur  at  the  interface.  The  main  idea  of  the  fast  Poisson  solver  on  an  irregular 
domain  is  to  extend  the  Poisson  equation  to  a  rectangular  domain  R  D  H.  This  procedure 
allows  us  to  use  fast  Poisson  solvers  on  a  fixed  Cartesian  grid  on  the  rectangular  domain,  for 
example,  the  Fishpack  [30]. 

The  fast  Poisson/Helmholtz  solvers  for  interior/exterior  problems  with  the  boundary  rep¬ 
resented  by  the  zero  level  set  of  a  two  dimensional  function  are  available  to  the  public  [18]. 
We  use  the  one  designed  for  interior  Helmholtz  equations.  The  package  also  provides  the 
orthogonal  projections  Xfc  of  the  interior  irregular  grid  points,  as  well  as  the  tangential  and 
normal  directions  of  the  boundary  dVL  at  those  projections. 

2.4  The  discrete  system  of  equations  in  matrix  vector  form 

Given  a  discrete  approximation  of  g(x,  y )  to  the  Laplacian  Au\gfi  along  the  boundary  at  those 
orthogonal  projections  X&,  we  can  use  the  fast  Poisson  solver  mentioned  above  to  solve  the 
two  Poisson  equations  in  (2.9)  to  get  ug(x,y).  We  denote  the  vector  of  the  discretized  values 
of  Uij  (from  the  interior  grid  points)  by  U;  and  the  vector  of  the  discretized  values  of  g(x,  y ) 
at  the  orthogonal  projections  of  the  interior  irregular  grid  points  by  G.  The  discrete  from 
(2.9)  can  be  written  as 


AXJ  +  BG  =  Fi 


(2.19) 
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for  some  vector  Fi  and  matrices  A  and  B1 .  It  requires  solving  two  Poisson  equations  on  the 
same  irregular  domain  11  with  different  Dirichlet  boundary  conditions  to  get  the  solution  U. 

Once  we  know  the  solution  U  of  the  augmented  system  (2.9)  given  G,  we  can  interpolate 
Utj  linearly  to  get  at  those  projections  X^,  1  <  k  <  JVj,.  The  interpolation  scheme  is 
crucial  to  the  success  of  our  algorithm  and  will  be  explained  in  detail  in  the  next  section. 
Therefore  we  can  write 


9U 


dn 


on 


CTJ  =  CA-1  (Fi  -  B G) 

Ck4_1Fi  -  CA~lBG, 


(2.20) 


where  C  is  a  sparse  matrix  determined  by  the  interpolation  scheme.  The  matrices  and  vectors 
are  only  used  for  theoretical  purposes  but  not  actually  constructed  in  our  algorithm.  We  need 
to  choose  such  a  vector  G  that  the  second  boundary  condition  |^|an  =  g2(x,y)  is  satisfied 
along  the  boundary  dil.  Therefore  we  have  the  second  matrix  vector  equation 

=  CA~1F1-CA~1BG  =  G2,  (2.21) 

an 

where  G2  is  the  vector  formed  from  the  boundary  condition  |)(|an  =  92{x ,  y)  at  X*.,  1  <  k  < 
Nfj.  Rewrite  the  equation  above  as 

EG  =  G2  -  CA~lF  1  =  F2,  (2.22) 


du 

dn 


where  E  =  —CA  lB  and  F2  =  G2  —  CA  1Fi.  If  we  put  the  two  matrix-vector  equations 
(2.19)  and  (2.22)  together  we  get 


'  A 

B  ' 

U  ' 

'  Fx  ' 

0 

E 

G 

.  F2  „ 

(2.23) 


Note  that  G  is  defined  only  on  a  set  of  points  X*.  while  U  is  defined  at  all  interior  grid  points. 
Let  TVjn  be  the  total  number  of  the  grid  points  in  the  interior  or  on  its  boundary  911,  then 
we  should  have  O(N^)  ~  0(lVjn).  The  Schur  complement  for  G  is 


EG  =  F2. 


(2.24) 


If  we  can  solve  for  the  system  above  to  get  G,  then  we  can  get  U  easily.  Because  the 
dimension  of  G  is  much  smaller  than  U,  we  expect  to  get  a  reasonably  fast  algorithm  for  the 
biharmonic  equation  on  irregular  domains  if  we  can  solve  (2.24)  efficiently. 

In  implementation,  we  use  the  GMRES  [25]  to  solve  (2.24).  The  GMRES  method  only 
requires  the  matrix  vector  multiplication.  We  explain  below  why  we  do  not  need  to  form  the 
matrix  E  explicitly. 

1Actually,  the  discrete  form  of  the  first  equation  in  (2.9)  can  be  written  as  TiV  =  F  +  EiG  for  some 
matrices  Ai  and  E\\  and  the  second  equation  can  be  written  as  TiU  =  V  +  F1G1.  Therefore  we  have 
A\  U  =  F  +  FiG  A1E1G1.  However,  the  matrices  .4 1 .  E \ ,  A,  and  B  are  used  only  for  theoretical  purposes 
and  never  formed  explicitly.  It  would  take  more  time  and  storage  to  compute  these  matrices  than  to  solve  the 
two  Poisson  equations  on  irregular  domain  directly. 


First  we  set  G  =  0  and  solve  the  two  Poisson  equations  in  (2.9),  or  (2.19)  in  the  discrete 
form,  to  get  U(0)  which  is  A~1F\  from  (2.19).  Note  that  the  residual  of  the  Schur  complement 
for  G  =  0  is 

R(  0)  =  F0-F2  =  -G2  +  £L4_1Fi 


=  -G2  +  CU(0)  =  -g2  + 


au(o) 


dn 


(2.25) 


an 


which  gives  the  right  hand  side  of  the  Schur  complement  system  with  an  opposite  sign.  The 
matrix-vector  multiplication  of  the  Schur  complement  system  given  G  is  obtained  from  the 
following  two  steps: 


Step  1:  Solve  the  two  Poisson  equations  in  (2.9),  or  (2.19)  in  the  discrete  form,  to  get  U(G). 
Step  2:  Interpolate  U(G)  to  get  |an-  Then  the  matrix  vector  multiplication  is 


EG 


aU(G) 


dn 


an 


<9U(0) 


On 


an 


This  is  because 


(2.26) 


EG  = 


from  the  equalities  E  = 


-CA^BG  =  - C  (A_1Fi  -  U) 


-CU(0)  +  CU(G)  = 


au(G) 


dn 


an 


au(o) 

dn 


—CA~1B,  AJJ  +  BG  =  Fi,  and  U(0)  =  A 


I  an 

_1Fi. 


Now  we  can  see  that  a  matrix  vector  multiplication  is  equivalent  to  solving  the  two  Poisson 
equations  in  (2.9),  or  (2.19)  in  the  discrete  form,  to  get  U,  and  using  an  interpolation  scheme 
to  get  Qj^\an  at  the  orthogonal  projections  of  the  interior  irregular  grid  points. 

While  our  approach  has  some  similarities  with  an  integral  equation  approach  to  find  a 
source  strength,  the  method  described  here  have  a  few  special  features:  (1)  no  Green  function 
is  needed  and  the  discussion  can  carry  over  to  three  dimensional  problems  directly;  (2)  no 
need  to  set  up  the  system  of  linear  equations  for  the  Schur  complement;  (3)  the  process 
itself  does  not  depend  on  the  boundary  condition  and  the  source  term.  The  efficiency  of  the 
algorithm  depends  on  the  number  of  iterations  of  the  GMRES  method. 


3  The  weighted  least  squares  interpolation  scheme 

The  interpolation  scheme  (2.20)  is  crucial  to  the  efficiency  (accuracy  and  the  number  of 
iterations  of  the  GMRES)  of  the  method.  Since  only  the  information  inside  the  domain  0  is 
useful,  the  interpolation  scheme  at  a  point  X  on  the  boundary  can  be  written  as 

E  7*j  Uij  ( | X  —  x* j  | ) ,  (3.27) 
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(3.28) 


where  da(r)  is  a  weighted  distance  function, 


dQ(r)  =  a5a/2(r) 


7^(1  +  cos(7rr/a)) 

0 


if  |r|  <  a 
if  Irl  >  a. 


Note  that  is  one  of  components  needed  in  the  matrix  vector  multiplication  EG.  We 

need  to  determine  the  parameter  a  and  the  interpolation  coefficients  7 ij.  Below  we  discuss 
how  to  determine  the  coefficients  7 jj.  Actually,  these  coefficients  are  different  from  point 
to  point  on  the  boundary.  So  they  should  really  be  labeled  as  But  for  simplicity  of 

notation,  we  will  concentrate  on  a  single  point  X  =  (A,  Y)  and  drop  out  the  subscript  X. 

Since  it  is  the  normal  derivative  that  we  are  trying  to  interpolate,  we  use  the  local  coor¬ 
dinates  at  the  boundary  point  (A,  Y), 


£  =  (x  —  A)  cos  0  +  (y  —  Y)  sin  6 , 
r]  =  — (x  —  A)  sin  9  +  (y  —  Y)  cos  9, 


(3.29) 


where  6  is  the  angle  between  the  x-axis  and  the  normal  direction  at  the  point  (A,  Y).  Under 
such  new  coordinates,  the  interface  can  be  parameterized  by  £  =  x{v)i  y  =  y.  Note  that 
x(0)  =  0,  and,  x'(0)  =  0  as  well,  assuming  that  the  boundary  is  smooth  enough  at  (A ,Y). 

We  use  an  un-determined  coefficients  method  to  determine  the  coefficients  7 ij.  Let  (£*,  rjj) 
be  the  £-77  coordinates  of  ( x%,yj ),  we  have  the  following  from  the  Taylor  expansion: 

u(xu  yj )  =  u(£i,  rjj)  =  u  +  +  uvyj  +  ^u^£2  +  ^um y]  +  u^yj  +  0(h3),  (3.30) 

where  for  simplicity,  we  use  the  same  notation  for  u  and  its  derivatives  in  the  original  and  the 
local  coordinates,  u,  u$,  •  ■  ■ ,  u^v  are  defined  at  (A,  Y)  in  the  original  coordinates,  or  (0, 0)  in 
the  local  coordinates.  Plugging  in  (3.30)  into  (3.27)  and  collecting  terms,  we  have 

~  ^2  7iJ  U(Xi ’  %')  da{\X  -  Xij\) 

0  (3.31) 

=  a\u  +  a2u^  +  a^Urj  +  04 +  asttw  +  clqu^  +  0(h 3  max  \jij  |) 
where  the  aj’s  are  given  by 


al  =  5^7ijda(|X-xy|) 

a4  —  ^2  |X  xy  ) 

=  'y  'J£k'Yijdg(\X-  xij  | ) 
ij 

“5  =  S  2r,^lijda^  ~  Xij^  (3-32) 

«3  =  'Y^rlklijda{  |X  -  Xij  |) 

l->3 

=  ^2  ^•kVk'yijdcti  |X  —  Xijl). 

i,j  i,j 


From  the  local  coordinate  transformation,  we  have  un  =  u%,  hence  we  can  set  up  the  linear 
system  of  equations  for  the  coefficients  7 y  as 


a-2  =  1, 
a5  =  0, 


ai  =  0 
04  =  0 
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a3  =  0, 

a6  =  0. 


(3.33) 


If  more  than  six  different  interior  grid  points  are  involved,  then  there  is  at  least  one  solution. 
Usually  we  choose  an  neighborhood  of  X  that  contains  more  than  six  different  interior  grid 
points  so  that  we  have  an  under-determined  system.  We  use  the  singular  value  decomposition 
(SVD)  to  find  the  least  squares  solution  with  the  least  1-2  norm.  In  this  way,  the  coefficients 
7 ij  have  roughly  the  same  magnitude  0(l/h),  and  7/)dQ.(|X  —  xy|)  is  roughly  a  decreasing 
function  of  the  distance  measured  from  X.  The  least  squares  interpolation  plays  an  important 
role  in  the  stability  of  the  algorithm.  In  practice,  only  a  hand  full  of  grid  points,  controlled  by 
the  parameter  a  and  the  normal  direction  at  the  boundary  point  ( X ,  Y),  are  involved.  Those 
grid  points  which  are  closer  to  ( X ,  Y)  have  more  influence  than  those  which  are  further  away. 

The  only  trade-off  of  our  weighted  least  square  interpolation  is  that  we  have  to  solve  an 
extra  under-determined  system  of  linear  equations.  The  larger  a  is,  the  more  computational 
cost  in  solving  (3.33).  However,  the  size  of  the  linear  system  is  small  and  the  coefficients  can 
be  pre-determined  before  the  GMRES  iteration.  We  will  see  that  the  extra  cost  in  dealing 
with  the  boundary  is  only  a  fraction  of  the  total  CPU  time  compared  to  the  Poisson  solvers. 
We  also  tried  the  third  order  interpolation  scheme  (10  equations)  and  the  numerical  results 
are  similar. 

3.1  Selecting  grid  points  for  the  interpolation  scheme 

If  the  interpolation  points  are  chosen  radially  from  the  center  of  the  interpolation  point 
X  until  enough  interior  grid  points  (6  ~  9  for  second  order  schemes,  10  ~  15  for  third 
order  scheme)  are  included,  the  method  works  fine  if  the  curvature  at  X  is  not  too  large. 
However,  if  the  curvature  is  large  and  there  are  fewer  grid  points  in  a  thin  tube  of  the  normal 
direction  compared  with  that  of  the  tangential  direction,  then  we  may  either  need  a  large 
circle  to  include  more  interior  points  in  the  tube,  or  we  may  not  have  a  good  accuracy  for 
the  interpolated  normal  derivative.  Either  case  will  affect  the  efficiency  of  the  numerical 
algorithm. 

Our  new  interpolation  scheme  is  to  select  the  grid  points  in  a  thin  cone.  The  vortex 
of  cone  is  the  interpolation  point  X  and  the  axis  is  the  normal  line  passing  though  X,  see 
Fig.  2  for  an  illustration,  and  see  [3]  for  more  details.  This  approach  works  better  because 
we  need  the  directional  derivative  of  u  along  the  normal  direction.  The  angle  of  the  cone  is 
a  parameter  such  that  tanfi  is  between  [1/10, 1/2]  in  our  choice. 

4  Numerical  examples 

We  have  done  a  number  of  numerical  experiments  which  confirm  the  expected  order  of  accu¬ 
racy.  All  the  computations  are  done  using  Sun  Ultra  10  workstations.  The  fast  Poisson  solver 
on  irregular  domain  used  is  from  the  IIM-pack  [18],  and  the  tolerance  for  the  convergence  is 
10~6  for  the  2-norm  of  the  residual  vector. 

Example  4.1  In  this  example  we  consider  a  biharmonic  equation  defined  on  a  circle  x2+y 2  = 
1/4  with  the  constructed  exact  solution 

u(x,  y)  =  x2  +  y2  +  ex,  (x,y)stl.  (4.34) 
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Figure  2:  Two  illustrative  examples  of  the  orthogonal  projection  (x*,y*)  of  irregular  grid 


points  on  the  boundary,  and  the  selected  grid  points  for  interpolation  scheme  at  the  projection 

toget^. 


The  forcing  term  f(x,  y)  is  obtained  by  applying  the  biharmonic  operator  to  the  exact  solution 
u(x,y ) 

f(x,y)  =  ex.  (4.35) 

The  normal  derivative  of  the  solution  on  the  boundary  is 

un(x,y)  =  8x2  +  4xex  +  8y2 ,  (x,y)&d£l,  (4.36) 

where  dTl  is  the  boundary  of  the  circle.  The  computation  domain  is  chosen  as  the  square 

[-1,1]  x  [-1,1]. 

Table  1  shows  the  results  of  a  grid  refinement  study,  where  the  first  column  is  the  number 

of  uniform  grid  points  in  the  x  and  y  directions.  The  maximum  error  is  defined  over  all  the 

interior  grid  points, 

||-En||oo  =  max  | u(xi,yj) -Uij\,  (4.37) 

h3,>P(hj)<  0 

where  Uij  is  the  computed  approximation  at  the  interior  grid  points  ( Xi,yj ).  The  third 
column  is  the  ratio  of  two  consecutive  errors  defined  as 

I  E  II 

|  -^n  1 1  oo 
hj2n  1 1  oo 

For  a  second  order  accurate  method,  the  ratio  should  approach  number  four  while  the  ratio 
should  approach  number  two  for  a  first  order  method.  We  can  see  clearly  second  order 
accuracy  of  our  method.  The  fourth  column  is  the  number  of  iterations.  We  see  that  only  a 
few  iterations  are  needed  and  the  number  is  independent  of  the  mesh  size.  The  fifth  column 
is  the  CPU  time  which  show  the  method  is  very  fast  considering  the  very  large  condition 
number  of  the  system  of  equations  from  a  direct  discretization  even  for  a  regular  domain 
such  as  rectangles.  We  will  use  the  same  notations  for  the  rest  of  examples  in  this  section. 


(4.38) 
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mesh  size 

\\En  |oo 

n 

No. 

CPU(s) 

64  x  64 

3.0003xl0”4 

4 

0.25 

128  x  128 

8.0076xl0”5 

3.74678 

5 

0.7709 

256  x  256 

1.9710xl0”5 

4.06270 

5 

2.3530 

512  x  512 

4.7654x10”° 

4.13604 

4 

10.375 

Table  1:  A  grid  refinement  analysis  for  example  (4.1). 


2  2 

Example  4.2  In  this  example,  the  boundary  is  a  skinny  ellipse  =  1.  The 

ential  equation  is: 

.  o  24y  12x  6.x3 

A  ^7/  —  -  —  -  —  - 

(1  +  x)5  (1  +  y)2  (1  +  y)4' 

We  use  the  Dirichlet  boundary  condition  which  is  determined  from  the  exact  solution 

u(x,y)  =  x3ln(l  +  y)  +  — y— ,  (x,y)eCl, 

1  +  x 


differ- 


(4.39) 


(4.40) 


and  the  normal  derivative  of  un(x,  y )  is  also  computed  from  the  exact  solution.  Note  that 
the  level  set  function  is  used  to  find  the  normal  and  tangential  directions.  The  computation 
domain  is  chosen  as  [—0.6,  0.6]  x  [—0.3, 0.3].  In  other  words,  we  can  use  a  rectangle  that  is 
just  large  enough  to  close  the  irregular  domain. 


Table  2  shows  the  results  of  a  grid  refinement  study.  We  can  see  that  the  method  still 
has  average  second  order  accuracy,  requires  only  a  few  iterations,  and  it  is  very  fast  in 
terms  of  the  CPU  time.  For  this  example,  the  curvature  is  quite  large  at  two  ends  of  the 
major  axes  of  the  ellipse.  The  new  interpolation  scheme  is  crucial  to  the  algorithm.  For 


mesh  size 

1  \En  |oo 

n 

No. 

CPU(s) 

64  x  32 

3.6576X  10”4 

8 

0.7310 

128  x  64 

9.5430X  10-5 

3.8327 

4 

1.0810 

256  x  128 

2.0889X  10~5 

4.5684 

5 

4.3360 

512  x  256 

4.9899x10”° 

4.1862 

6 

21.000 

Table  2:  A  grid  refinement  analysis  for  example  (4.2). 


the  biharmonic  equations,  the  accuracy  of  solution  depends  on  the  relative  position  of  the 
irregular  boundary  and  the  underlying  Cartesian  grid,  see  [17].  It  is  more  reasonable  to  do 
the  linear  regression  analysis  to  find  the  convergence  order  which  is  shown  in  Fig.  3.  The 
average  order  of  convergence  is  2.0240. 

Example  4.3  In  this  example  we  consider  a  biharmonic  equation  defined  on  a  non-convex, 
non-concave  region  defined  by 

f  x=  (0.6  +  0.25  sin(5#))  cos(0), 

<  0e  0,27t).  (4.41) 

!  y=  (0.6  +  0.25  sin(50))  sin(0), 
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Figure  3:  The  linear  regression  analysis  on  the  errors  for  example  (4.2). 


mesh  size 

hr  ii 

1 1  \  \  oo 

n 

No. 

CPU(s) 

64  x  64 

8.14721xl0”4 

9 

2.2130 

128  x  128 

2.08149x  10~4 

3.91412 

7 

3.7360 

256  x  256 

8.06897xl0-5 

2.57962 

9 

13.520 

512  x  512 

1.83151xl0-5 

4.40563 

9 

52.034 

Table  3:  A  grid  refinement  analysis  for  example  (4.3) 


The  differential  equation  is: 

A  2u  =  0.  (4.42) 

The  essential  boundary  condition  u\g q  and  un\gn  are  determined  from  the  constructed  exact 
solution 

u(x,  y)  =  x2  +  y2  +  excos(y),  (4.43) 

in  connection  with  the  level  set  function. 

Table  3  shows  the  results  of  a  grid  refinement  study.  Fig.  4  shows  the  mesh  plot  of  the 
exact  solution.  Fig.  5  shows  the  linear  regression  analysis.  The  average  convergence  order  is 
1.9300. 

4.1  Algorithm  efficiency  analysis 

We  have  seen  from  previous  examples  that  the  number  of  iterations  are  fairly  small.  For 
most  of  the  examples,  it  is  between  4  ~  10.  Another  concern  about  the  algorithm  is  how 
much  overhead  cost  is  needed  for  dealing  with  the  interior  irregular  points.  The  cost  includes 
indexing  the  irregular  grid  points,  finding  the  projections,  and  solving  a  linear  system  of 
equations  to  find  the  coefficients  of  the  interpolation  scheme  at  each  orthogonal  projection 
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Figure  4:  The  computed  solution  for  example  (4.3).  The  difference  between  the  exact  and 
computed  solutions  is  hardly  visible. 


-8  -7  -6  -5  -4  -3  -2 


logl  0(h) 


Figure  5:  The  linear  regression  analysis  on  the  errors  for  example  (4.3). 
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on  the  boundary  of  the  interior  irregular  points.  Certainly  the  CPU  time  depends  on  the 
geometry.  In  Table  4,  we  show  the  CPU  time  spent  in  dealing  with  the  boundary  and  that 
in  solving  the  entire  problem.  We  can  see  that  the  overhead  time  is  only  a  fraction  of  the 
total  CPU  time  especially  as  the  mesh  gets  finer. 


mesh  size 

Example  (4.1) 

Example  (4.3) 

tov 

tsoive 

tov 

tsoive 

64  x  64 

0.0499 

0.2510 

0.0699 

2.2130 

128  x  128 

0.1100 

0.7749 

0.2710 

3.7360 

256  x  256 

0.1210 

2.3530 

0.5410 

13.520 

512  x  512 

0.3999 

10.3750 

0.6610 

52.035 

Table  4:  The  CPU  time  in  dealing  with  interior  irregular  grid  points  (toy)  and  the  CPU  time 
for  the  linear  solver  (tsoive)  for  Example  (4.1)  and  (4.3). 


5  An  applications  of  the  biharmonic  solver  to  the  incompress¬ 
ible  Stokes  flow  on  an  irregular  domain 

Consider  the  incompressible  Stokes  flow  in  two  dimensional  spaces 

pAu  =  px  —  F\ ,  x  G  P, 

pAv  =  py  —  F‘2 ,  x  G  0,  (5.44) 

V  •  u  =  0, 

where  p  is  the  fluid  viscosity,  p  is  the  pressure,  u  =  (u,  v)  is  the  velocity,  and  F  =  (T\,  F2)  is 
a  forcing  term.  Equations  (5.44)  are  supplemented  by  the  “no-slip”  boundary  condition 

u(x)|an  =  0.  (5.45) 

The  vorticity  function,  which  is  a  scalar  in  two-dimensional  case,  is  defined  by 

ca(x)  =  (V  x  u)  •  k  =  —uy  +  vx,  x  6  P,  (5.46) 

where  k  is  the  unit  vector  in  the  ^-direction.  The  first  two  equation  of  (5.47)  can  be  written 
in  vector  form 

pAu  =  Vp  —  F,  xeQ.  (5-47) 

Taking  curl  of  the  above  equation  we  get 

—pAu:  =  (V  x  F)  •  k,  xell,  (5.48) 

The  velocity  field  u  is  obtained  by  using  (5.46)  in  conjunction  with  (5.45).  This  can  be  done 
via  the  stream-function  formulation  as  follows.  Due  to  V  •  u  =  0,  there  is  a  scalar  function 
V’(x)  such  that 

u(x)  =  V±V’=  (“§|>^)>  xGa  (5-49) 
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Therefore  from  (5.46),  we  get 


Ai/j  =  uj.  (5.50) 

Thus,  equations  (5.48)  and  (5.50),  along  with  the  relation  (5.49),  serve  as  the  vorticity  stream- 
function  formulation  of  the  problem. 

We  note  that  the  “no-slip”  condition  (5.45)  and  the  relation  (5.49)  implies 


V’Wlan  =  J^(x) 


=  0, 


an 


(5.51) 


where  n  is  the  unit  normal  vector  of  the  boundary  <917  pointing  outward.  Note  that  (5.51) 
follows  since  if}  is  a  constant  along  the  boundary  and  is  only  determined  up  to  an  additive 
constant. 

Plugging  in  (5.50)  into  (5.48),  and  using  the  relation  (5.51),  we  get 


— /iA2^  =  (V  x  F)  •  k,  xe!l, 
V>(x)lan  =  °> 

^n(x)|an  =  0. 


(5.52) 


This  is  a  well  defined  biharmonic  equation  with  essential  boundary  condition.  The  advantage 
of  using  the  formulation  to  solve  for  the  velocity  field  u  is  that  it  eliminates  the  difficulty  of 
dealing  with  the  improper  partition  of  boundary  conditions  of  the  vorticity-stream  function 
formulation.  As  we  stated  earlier,  we  assume  that  the  boundary  of  the  domain  <911  is  piecewise 
smooth. 

We  use  the  algorithm  developed  in  previous  sections  to  solve  (5.52)  numerically.  Once  V’ 
is  obtained,  the  standard  central  difference  formula  can  be  used  to  interpolate  (5.49)  to  get  u 
and  v  at  regular  grid  points.  For  irregular  interior  grid  points,  we  use  a  similar  least  squares 
interpolation  technique  developed  in  Section  3  to  recover  u  and  v. 

We  have  tested  our  algorithm  for  the  Stokes  flow  with  Reynolds  number  being  400.  The 
stream  function  is 

y)  =  -  sin2(X  +  V  tt)  -  -.  (5.53) 

7T  4  7T 

The  fluid  is  confined  in  a  circular  region 

&  =  {(x,y)\  x2 +  y2  =  2}  (5.54) 


and  the  computation  region  is  chosen  to  be 
condition 


ip(x,y)  =  0, 


n(x,y )  =  0, 


[—1.6, 1.6]  x  [—1.6, 1.6].  The  no-slip  boundary 


(x,y)  €  <911, 
( x ,  y  )  e  <911. 


(5.55) 


are  satisfied. 

In  Table  5,  we  show  the  results  of  a  grid  refinement  study  for  the  velocity  field  u  and  v. 
We  can  see  clearly  second  order  accuracy  for  the  solution  ( u ,  v)  in  the  infinity  norm.  We  also 
tested  other  examples  and  got  similar  results. 
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mesh  size 

1 1 E  (r)  n  1 1  oo 

r 

||  ||oo 

r 

64  x  64 

1.07739x  10-2 

1.08006  x  10“2 

128  x  128 

2.54015x  10-3 

4.2415 

2.56114  x  10“3 

4.2171 

256  x  256 

4.49293x  10-4 

5.6570 

4.49161  x  10~4 

5.7037 

512  x  512 

1.09871xl0-4 

4.0892 

1.09842  x  10“4 

4.0907 

Table  5:  A  grid  refinement  analysis  of  u  =  (u,v)  for  example  (5.53)-(5.55). 


6  Summary  and  acknowledgment 

In  this  paper,  a  fast  iterative  method  has  been  developed  for  biharmonic  equations  defined  on 
irregular  domains  with  essential  boundary  conditions.  The  method  takes  advantage  of  the  fast 
Poisson  solver  on  irregular  domains  available  to  the  public  and  uses  an  augmented  approach  so 
that  the  biharmonic  equation  is  decomposed  to  two  Poisson  equations.  The  GMRES  method 
is  used  to  solve  the  augmented  variable  which  is  the  Laplacian  of  the  solution  along  the 
boundary.  An  interpolation  scheme  for  approximating  the  normal  derivative  of  the  solution 
is  crucial  to  the  success  of  the  algorithm.  Numerical  examples  show  the  efficiency  of  the  second 
order  accurate  method.  Only  a  few  iterations  are  needed  to  solve  the  augmented  variable 
and  the  entire  problem.  An  application  of  the  fast  biharmonic  solver  to  the  incompressible 
Stokes  flow  is  also  presented. 
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